################################################################################
########################09. CSES Only##################################
################################################################################

pols <- readRDS('pols_bjps')

#####Load plm package#####
#install.packages('plm')
library(plm)
#install.packages('pcse')
library(pcse)
#install.packages('lmtest')
library(lmtest)
#install.packages('dplyr')
library(dplyr)


###Fixed Effects Model (1)###

#####CSES only#####

cses <- plm(scale_econsdw ~ scale_dispgini + 
              scale_ggini_inc_cs + 
              scale_socioeconomicsal +
              scale_dispgini : scale_ggini_inc_cs + 
              scale_dispgini*scale_socioeconomicsal +
              scale_ggini_inc_cs: scale_socioeconomicsal +
              scale_ggini_inc_cs:scale_socioeconomicsal : scale_dispgini +
              scale_galtansdw +
              scale_enep +
              scale_unemp,
            data = pols, 
            na.action = na.omit, 
            model = "within")

print(coeftest(cses, vcov = function(x) vcovBK(x, cluster="group")))

print(summary(cses) )#Model Fit Stats



#####NOT PRESENTED IN TEXT#####
###Multi-Level Model (2)###

#install.packages('nlme')
library(nlme)

####Only CSES Data####

csese <- lme(scale_econsdw ~  scale_ggini_inc_cs +
               scale_socioeconomicsal +
               scale_ggini_inc_cs*scale_socioeconomicsal +
               scale_dispgini+
               scale_enep +
               scale_galtansdw +
               scale_mag+
               scale_unemp +
               scale_ggini_inc_cs*scale_dispgini +
               scale_dispgini * scale_socioeconomicsal +
               scale_ggini_inc_cs * scale_dispgini * scale_socioeconomicsal,
             random = ~ 1 | eastwest / country.name , 
             data = pols,
             na.action = na.omit,
             method = 'ML')

summary(csese)

#Model Fit Statistics#

stargazer::stargazer(csese)


library(MuMIn)

print(r.squaredGLMM(csese))

print(summary(csese)$logLik)
print(summary(csese)$AIC)
print(summary(csese)$BIC)